Modélisation du double pendule

Lagrangien Libre

Le lagrangien pour deux particules de masse $m_i$ de coordonnées $(x_i, y_i)$ dans un champs gravitationnels constant, (l'ordonnée pointant vers le bas) vaut $$L:=m_1\frac{\dot{x_1}^2+\dot{y_1}^2}{2}+m_1gy_1+m_2\frac{\dot{x_2}^2+\dot{y_2}^2}{2}+m_2gy_2.$$

Contrainte cinématique

Dans le cas du double pendule de longueur $l_i$ on introduit les trois angles entre la verticale et le pendule et on obtient $$ \begin{cases} x_1=r_1\sin(\theta_1)\\ y_1=r_1\cos(\theta_1) \end{cases},\qquad \begin{cases} x_2=x_1+r_2\sin(\theta_2)\\ y_2=y_1+r_2\cos(\theta_2) \end{cases},\qquad $$

Lagrangien nouvelles coordonnées

Le lagrangien s'écrit alors $$ L:=m_1\frac{(r_1 \dot{\theta_1})^2}{2}+m_1gr_1\cos(\theta_1)+m_2\left(\frac{(r_1 \dot{\theta_1})^2}{2}+\frac{(r_2 \dot{\theta_2})^2}{2}+r_1\dot{\theta_1}r_2\dot{\theta_2}\cos(\theta_1-\theta_2)\right)+m_2g\left(r_1\cos(\theta_1)+r_2\cos(\theta_2)\right) $$

Equations d'Euler-Lagrange

On peut maintenant écrire les équations d'Euler-Lagrange $$ \begin{cases} \frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L}{\partial \dot{\theta_1}}=\frac{\partial L}{\partial\theta_1}\\ \frac{\mathrm{d}}{\mathrm{d}t}\frac{\partial L}{\partial \dot{\theta_2}}=\frac{\partial L}{\partial\theta_2} \end{cases} $$

Or on peut calculer $$ \begin{cases} \frac{\partial L}{\partial \dot{\theta}_1}=m_1 r_1^2\dot{\theta}_1+m_2 r_1^2\dot{\theta}_1+m_2r_1r_2 \dot{\theta}_2\cos(\theta_1-\theta_2)\\ \frac{\partial L}{\partial \dot{\theta}_2}=m_2 r_2^2\dot{\theta}_2+m_2r_1r_2 \dot{\theta}_1\cos(\theta_1-\theta_2)\\ \frac{\partial L}{\partial \theta_1}=-m_1gr_1\sin(\theta_1)-m_2r_1r_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_2gr_1\sin(\theta_1),\\ \frac{\partial L}{\partial \theta_2}=m_2r_1r_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_2gr_2\sin(\theta_2), \end{cases} $$

Ce qui donne donc les équations d'Euler-Lagrange suivantes $$ \begin{cases} m_1 r_1^2\ddot{\theta}_1+m_2 r_1^2\ddot{\theta}_1+m_2r_1r_2 \ddot{\theta}_2\cos(\theta_1-\theta_2)+m_2r_1r_2 (\dot{\theta}_2)^2\sin(\theta_1-\theta_2)-m_2r_1r_2\dot{\theta}_1 \dot{\theta}_2\sin(\theta_1-\theta_2)=-m_1gr_1\sin(\theta_1)-m_2r_1r_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_2gr_1\sin(\theta_1)\\ m_2 r_2^2\ddot{\theta}_2+m_2r_1r_2 \ddot{\theta}_1\cos(\theta_1-\theta_2)-m_2r_1r_2 (\dot{\theta}_1)^2\sin(\theta_1-\theta_2)+m_2r_1r_2 \dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)=m_2r_1r_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_2gr_2\sin(\theta_2) \end{cases} $$

Et après rapide simplification $$ \begin{cases} (m_1+m_2) r_1^2\ddot{\theta}_1+m_2r_1r_2 \ddot{\theta}_2\cos(\theta_1-\theta_2)=-(m_1+m_2)gr_1\sin(\theta_1)-m_2r_1r_2 (\dot{\theta}_2)^2\sin(\theta_1-\theta_2)\\ m_2 r_2^2\ddot{\theta}_2+m_2r_1r_2 \ddot{\theta}_1\cos(\theta_1-\theta_2)=m_2r_1r_2 (\dot{\theta}_1)^2\sin(\theta_1-\theta_2)-m_2gr_2\sin(\theta_2) \end{cases} $$

Mise sous forme définie

On voit que les termes de plus haut de degrés apparaissent linéairement. La matrice étant $$ M(\theta_1,\theta_2):=\begin{pmatrix} (m_1+m_2)r^2_1 & m_2r_1r_2\cos(\theta_1-\theta_2)\\ m_2r_1r_2\cos(\theta_1-\theta_2) & m_2r_2^2 \end{pmatrix} $$

Or le déterminant de celle-ci est $$d:=\mathrm{det}M(\theta_1,\theta_2)=m_2 (r_1r_2)^2(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))>0.$$

On arrive alors à $$ M^{-1}(\theta_1, \theta_2)= \frac{1}{m_2 (r_1r_2)^2(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))} \begin{pmatrix} m_2 r_2^2 & -m_2r_1r_2\cos(\theta_1-\theta_2)\\ -m_2r_1r_2\cos(\theta_1-\theta_2) & (m_1+m_2)r_1^2 \end{pmatrix} $$

On peut alors l'inverser pour obtenir le système \begin{equation} \begin{pmatrix} \ddot{\theta}_1\ \ddot{\theta}_2 \end{pmatrix} = M^{-1}(\theta_1, \theta_2) \begin{pmatrix} -(m_1+m_2)gr_1\sin(\theta_1)-m_2r_1r_2 (\dot{\theta}_2)^2\sin(\theta_1-\theta_2)\ m_2r_1r_2 (\dot{\theta}_1)^2\sin(\theta_1-\theta_2)-m_2gr_2\sin(\theta_2) \end{pmatrix} \end{equation}

Formulation Hamiltonienne

Changement de variable

On pose $$ \begin{cases} p_1:=\frac{\partial L}{\partial \dot{\theta_1}}=m_1 r_1^2\dot{\theta}_1+m_2 r_1^2\dot{\theta}_1+m_2r_1r_2 \dot{\theta}_2\cos(\theta_1-\theta_2)\\ p_2:=\frac{\partial L}{\partial \dot{\theta_2}}=m_2 r_2^2\dot{\theta}_2+m_2r_1r_2 \dot{\theta}_1\cos(\theta_1-\theta_2) \end{cases} $$

On voit que la relation s'écrit aussi $$ \begin{pmatrix} p_1 \\ p_2 \end{pmatrix} = M(\theta_1, \theta_2) \begin{pmatrix} \dot{\theta}_1 \\ \dot{\theta}_2 \end{pmatrix} $$

et donc $$ \begin{cases} \dot{\theta}_1=\frac{m_2r_2^2p_1-m_2r_1r_2\cos(\theta_1-\theta_2)p_2}{m_2 (r_1r_2)^2(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))}\\ \dot{\theta}_2=\frac{(m_1+m_2)r_1^2p_2-m_2r_1r_2\cos(\theta_1-\theta_2)p_1}{m_2 (r_1r_2)^2(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))} \end{cases} $$

Hamiltonien

Mais on sait alors que si l'on construit la fonction $$H:=p_1\dot{\theta}_1+p_2\dot{\theta}_2-L(\theta_1,\theta_2,\dot{\theta}_1,\dot{\theta}_2)= (m_1+m_2)\frac{(r_1 \dot{\theta_1})^2}{2}+m_2\frac{(r_2 \dot{\theta_2})^2}{2}+m_2r_1r_2\dot{\theta_1}\dot{\theta_2}\cos(\theta_1-\theta_2)-(m_1+m_2)gr_1\cos(\theta_1)-m_2 g r_2\cos(\theta_2) ,$$ On a alors $$ \begin{cases} \dot{\theta}_1=\frac{\partial H}{\partial p_1}\\ \dot{\theta}_2=\frac{\partial H}{\partial p_2}\\ \dot{p}_1=-\frac{\partial H}{\partial \theta_1}\\ \dot{p}_2=-\frac{\partial H}{\partial \theta_2} \end{cases} $$

Or on a $$\dot{\theta}_1^2=\frac{1}{d^2}\left(m_2^2r_2^4p_1^2+m_2^2r_1^2r_2^2\cos^2(\theta_1-\theta_2)p_2^2-2m_2^2 r_1 r_2^3\cos(\theta_1-\theta_2)p_1 p_2\right).$$ $$\dot{\theta}_2^2=\frac{1}{d^2}\left((m_1+m_2)^2r_1^4p_2^2+m_2^2r_1^2r_2^2\cos^2(\theta_1-\theta_2)p_1^2-2m_2(m_1+m_2) r_1^3 r_2\cos(\theta_1-\theta_2)p_1 p_2\right).$$ $$\dot{\theta}_1\dot{\theta}_2=\frac{1}{d^2}\left((m_2(m_1+m_2)r_1^2r_2^2+m_2^2r_1^2r_2^2\cos^2(\theta_1-\theta_2))p_1p_2-m_2^2r_1r_2^3\cos(\theta_1-\theta_2)p_1^2-m_2(m_1+m_2)r_1^3r_2\cos(\theta_1-\theta_2)p_2^2 \right).$$

On peut donc écrire $$p_1\dot{\theta}_1+p_2\dot{\theta}_2-L(\theta_1,\theta_2,\dot{\theta}_1,\dot{\theta}_2)= (m_1+m_2)\frac{(r_1 \dot{\theta_1})^2}{2}+m_2\frac{(r_2 \dot{\theta_2})^2}{2}+m_2r_1r_2\dot{\theta_1}\dot{\theta_2}\cos(\theta_1-\theta_2)=\frac{1}{d^2}\Big(\\ p_1^2\left(\frac{(m_1+m_2)m_2^2}{2}r_1^2r_2^4+\frac{m_2^3}{2}r_1^2r_2^4\cos^2(\theta_1-\theta_2) -m_2^3r_1^2r_2^4\cos^2(\theta_1-\theta_2)\right)\\ +p_2^2\left(\frac{(m_1+m_2)m_2^2}{2}r_1^4r_2^2\cos^2(\theta_1-\theta_2)+\frac{m_2(m_1+m_2)^2}{2}r_1^4r_2^2 -m_2^2(m_1+m_2)r_1^4r_2^2\cos^2(\theta_1-\theta_2)\right)\\ p_1p_2\left(m_2^2r_1^3r_2^3\cos(\theta_1-\theta_2)(m_1+m_2(1+\cos^2(\theta_1-\theta_2))) -2(m_1+m_2)m_2^2r_1^3r_2^3\cos(\theta_1-\theta_2)\right) \Big)\\ =\frac{1}{d^2}\Big(p_1^2m_2^2r_1^2r_2^4\frac{m_1+m_2(1-\cos^2(\theta_1-\theta_2))}{2}\\ p_2^2(m_1+m_2)m_2r_1^4r_2^2\frac{m_1+m_2(1-\cos^2(\theta_1-\theta_2))}{2}\\ -p_1p_2m_2^2r_1^3r_2^3\cos(\theta_1-\theta_2)(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))\Big) $$

On peut alors écrire le Hamiltonien $$ H=\frac{m_2r_2^2p_1^2+(m_1+m_2)r_1^2p_2^2+m_2r_1r_2\cos(\theta_1-\theta_2)p_1p_2}{m_2r_1^2r_2^2(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))}-(m_1+m_2)gr_1\cos(\theta_1)-m_2gr_2\cos(\theta_2).$$

Schéma Numérique

On peut commencer par constater que l'on n'a pas besoin de $H$ pour obtenir une discrétisation car avoir accès aux valeurs successives de $\theta_1$ et $\theta_2$ permet d'obtenir les valeurs approchées de $\dot{\theta}_1$ et $\dot{\theta}_2$ or $$ \begin{cases} \dot{\theta}_1=\frac{m_2r_2^2 p_1-m_2r_1r_2\cos(\theta_1-\theta_2)p_2}{m_2 r_1^2r_2^2(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))}\\ \dot{\theta}_2=\frac{(m_1+m_2)r_1^2 p_2-m_2r_1r_2\cos(\theta_1-\theta_2)p_1}{m_2 r_1^2r_2^2(m_1+m_2(1-\cos^2(\theta_1-\theta_2)))}\\ \dot{p}_1=-m_1gr_1\sin(\theta_1)-m_2r_1r_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_2gr_1\sin(\theta_1),\\ \dot{p}_2=m_2r_1r_2\dot{\theta}_1\dot{\theta}_2\sin(\theta_1-\theta_2)-m_2gr_2\sin(\theta_2), \end{cases} $$

import numpy as np
import matplotlib.pyplot as plt
import  matplotlib.animation as animation
%matplotlib inline
plt.rc('figure', figsize=(12,9))

plt.rcParams['animation.writer'] = "avconv"

# Constantes du système
g, m1, m2, r1, r2 = 9.8, 2, 1, 1, 2

# flux pour les angles
def fluxth(th1, th2, p1, p2):
    c1 = (m2*r2**2*p1-m2*r1*r2*np.cos(th1-th2)*p2)/(m2*r1**2*r2**2*(m1+m2*(np.sin(th1-th2))**2))
    c2 = ((m1+m2)*r1**2*p2-m2*r1*r2*np.cos(th1-th2)*p1)/(m2*r1**2*r2**2*(m1+m2*(np.sin(th1-th2))**2))
    return c1, c2

# flux pour les quantités de mouvements
def fluxp(th1, th2, dth1, dth2):
    c1 = -(m1+m2)*g*r1*np.sin(th1)-m2*r1*r2*dth1*dth2*np.sin(th1-th2)
    c2 = m2*r1*r2*dth1*dth2*np.sin(th1-th2)-m2*g*r2*np.sin(th2)
    return c1, c2    

# Schéma d'Euler explicite
def euler(th10, th20, p10, p20, dt=0.01, T=10):
    N = int(T/dt)
    temps = np.linspace(0, T, N+1)
    sol = np.zeros((4,N+1))
    sol[0,0] = th10
    sol[1,0] = th20
    sol[2,0] = p10
    sol[3,0] = p20
    for i in range(N):
        a, b = fluxth(sol[0,i], sol[1,i], sol[2,i], sol[3,i])
        sol[0,i+1] = sol[0,i] + dt*a
        sol[1,i+1] = sol[1,i] + dt*b
        c, d = fluxp(sol[0,i], sol[1,i], a, b)
        sol[2,i+1] = sol[2,i] + dt*c
        sol[3,i+1] = sol[3,i] + dt*d
    return temps, sol

# Affichage
def visualisation(temps, sol):
    fig, ax = plt.subplots()
    x1 = r1*np.sin(sol[0,:])
    y1 = -r1*np.cos(sol[0,:])
    x2 = x1 + r2*np.sin(sol[1,:])
    y2 = y1 - r2*np.cos(sol[1,:])
    ax.set_xlim(-(r1+r2), r1+r2)
    ax.set_ylim(-(r1+r2), r1+r2)
    line, = plt.plot([0,x1[0], x2[0]], [0, y1[0], y2[0]], linewidth=2, marker="o", color="blue")
    def maj(i):
        line.set_data([0,x1[i], x2[i]], [0, y1[i], y2[i]])
        return line,
    
    ani = animation.FuncAnimation(fig, maj, frames=range(np.size(x1)), blit=True, interval=50)
    return ani

# application

temps, sol = euler(np.pi/2, np.pi/4, 1.5, 0, 0.05, 20)
ani = visualisation(temps, sol)

from IPython.display import HTML

HTML(ani.to_html5_video())